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Abstract: We propose a lattice model to study the dynamics of a driven 
interface in a medium with random pinning forces. The critical exponents 
characterizing the depinning transition are determined numerically in 1+1 
and 2+1 dimensions. Our findings are compared with recent numerical and 
analytical results for a Langevin equation with quenched noise, which is ex- 
pected to be in the same universality class. 
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I. Introduction 

In recent years there has been considerable interest in the scaling behaviour of 
interfaces and directed lines broadened by randomness [1] and in the critical phenom- 
ena which occur at the onset of steady-state motion [2,3]. The driven viscous motion 
of an interface in a medium with random pinning forces combines aspects of both 
of these two fields. An example is the motion of a domain wall in the random field 
Ising model [4,5]. The relaxation of metastable domains in a diluted antiferromagnet 
with applied external magnetic field provides a possible experimental realization [6] . 
Due to random field pinning the domain walls become rough and if the driving field 
is sufficiently weak compared to the random field strength, the domain walls become 
stuck and the domain state is frozen. Random field pinning can also be important for 
fluid displacement experiments in porous media [7], where an anomalous roughness 
in comparison with the KPZ-equation [1] was observed. Kessler, Levine and Tu [8] 
(KLT) and recently He, Kahanda and Wong [7] suggested that the increased rough- 
ness is due to quenched random capillary forces on the interface. (Another interesting 
proposal [9] is to assume that the noise fluctuates in time but follows a power-law 
distribution. This possibility will not be considered here.) Very similar problems 
arise in vortex pinning in type-II superconductors [10] and in sliding charge- density 
waves [2,11]. 

The simplest continuum description for the driven motion of an interface in a 
random medium is given by the following Langevin equation for a D-dimensional 
interface profile z(x, t) [12,13] 

\^ = ^V 2 z + F + r ] (x,z) (1) 

where A and 7 are the inverse mobility and the stiffness constant, respectively, and 
F is a uniform driving force. The random force 77 (x, z) is Gaussian distributed with 
{rfj = and correlations (r/(xo, zo)??(xo + x, z + z)) = S D (x.)A(z). For the random 
field case the correlator A(z) is a monotonically decreasing function of z for z > and 
decays rapidly to zero over a finite distance. The difference to the Edwards- Wilkinson 
(EW) equation [14] is that the noise depends on the position of the interface, which 
makes the problem highly nonlinear. However, the noise that acts on a sufficiently 
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fast moving interface is fluctuating in time like in the EW equation. Therefore a 
crossover to a EW regime is expected as the velocity increases to large values. Here 
we are interested in the critical behaviour at the onset of motion where the external 
driving force is just able to overcome the pinning forces. The velocity v of the interface 
scales with the reduced force f = (F — F c )/F c as v ~ f e where F c is the threshold 
force [2,15]. It can be shown [15] that there is a correlation length £ which diverges 
with F — > F c as £ ~ f~ u . Thus, at F = F c , the roughness (or interface width) 



where L is the system size, ( is the roughness exponent, z is the dynamical exponent 
and ^ is a scaling function with *ff(y) ~ y^ (f3 = (/z) for y <C 1 but becomes 
constant for y ^> 1. In eq.(2) the overbar denotes the spatial average over x and the 
angular brackets mean the configurational average. It is known [12,15] that above 
four interface dimensions D the interface is flat, i.e. £ = 0. For D = 4 — e, e < 1 
the critical exponents were calculated by a functional renormalization group scheme 
to first order in e: C - e/3, z = (/(3 ~ 2 - 2e/9, i/ ~ 2/(6 - e) and 9 = v{z - Q 
[15]. The equation of motion (1) was simulated by KLT in 1+1 dimensions for a 
wide range of velocities [8]. For v — > they suggest that C — > 1 which would be in 
agreement with an Imry-Ma argument [4]. Parisi [17] proposed that the argument 
z of the random force rj in eq.(l) can be replaced by a constant from which follows 
C = (4 - D)/2, (3 = (4 - D)/4 and = 1 for D < 4. He numerically determined 
(3 in D = 1,2 and 3 from (1) with results consistent with his analytical arguments. 
Recently, Tang [18] used the Runge-Kutta scheme to solve eq.(l) and found in 1+1 
dimensions C = 1-25 ± 0.10, (3 = 0.81 ± 0.02, 9 = 0.4 ± 0.05 and v = 1.1 ± 0.1, in 
disagreement with both KLT [8] and Parisi [17]. 



Here we introduce a lattice model of probabilistic cellular automata [19] which 
is expected to be in the same universality class as eq.(l). In 1+1 dimensions the 
results are ((D = 1) = 1.25 ± 0.01, (3{D = 1) = 0.88 ± 0.02, and in 2+1 dimensions 
C(D = 2) = 0.75 ± 0.02, (3(D = 2) = 0.475 ± 0.015 and 9(D = 2) = 0.65 ± 0.05. 




(2) 



scales as [16] 




(3) 
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These values support the results of the e-expansion [15] because the deviations for 
4 — D = e = 2 are much smaller than for e = 3. 

II. The model 

In order to simplify the notation the model is explained for the case of 1+1 
dimensions. The generalization to higher dimensions is straightforward. Consider 
a square lattice where each cell (i, h) is assigned a random pinning force 77^ which 
takes the value 1 with probability p and -1 with probability q = 1 — p. By excluding 
overhangs the interface is specified by a set of integer column heights hi(t), i = 
1, 2...L. At t = all columns have the same height hi(t = 0) = 0. During the motion 
for a given time t the value 

Vi = h i+1 (t) + hi-!(t) - 2hi(t) + grji,h (4) 

is determined for all i, where g is a parameter, 77^ = ±1 and periodic boundary 
conditions are used. The interface configuration is then updated simultaneously for 
all i: 

hi(t+l) = hi(t) + l if Vi>0 

hi(t+ 1) = hi(t) otherwise. (5) 

The parameter g measures the strength of the random force compared to the elastic 
force and the difference p — q determines the driving force. The growth rule specified 
by eqs.(4) and (5) can be motivated by the continuum equation (1) because Vi is the 
sum of the discretized Laplacian and a random force. It remains however an open 
question whether the use of a non-Gaussian noise and the discretization of x, z and t 
changes the universality class on accessible length scales. The relation of this model 
to the mechanism introduced earlier [20] for pinning by directed percolation will be 
discussed below. 

III. Numerical results 

a) 1+1 dimensions 
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Since we are interested in the scaling behaviour at the depinning transition one 
has to find first the critical value p c (g) of the concentration of plus cells where the 
infinite system gets pinned. For g = 1 we estimate p c (g = 1) ~ 0.8004, while 
p c (g = 2) ~ 0.8748. To determine the roughness exponent C, different system sizes 
with 8 < L < 1024 were simulated at p = p c until they get pinned. The values 
for w 2 (L) were averaged over several thousand independent runs n, depending on L 
( n\fL ~ 8* 10 4 ). Fig.l shows w 2 (L) in a double logarithmic plot. For 32 < L < 512 
the data are well fitted by a straight line from which we get ((D = 1) = 1.25 ± 0.01. 
Error bars indicate only statistical uncertainities. The fluctuations are somewhat 
larger for g = 2 than for g = 1. For bigger L the simulations become very time 
consuming because it takes longer for a larger system to get pinned at p = p c . In 
addition, when one increases L, one has to know the threshold p c more accurately 
to assure £ 3> L. To see the deviations from the scaling behaviour we define an 
effective roughness exponent C(-^) = l°g[iu(2L)/-u;(L)]/ log2. In fig. 2 ((L) for g = 1 
and p = p c is compared with ((L) for g = 2 and p — 0.8744 < p c . For L = 256 and 
512, C(L) is significantly smaller for p < p c than for p = p c . 

Besides the possibility to determine the dynamical roughness exponent (3 from the 
width w(t) it is also possible to consider the scaling behaviour of H(t) = {hi) ~ t 13 . 
From the latter we get (3(D = 1) = 0.86 ± 0.02 for L = 16384, while from the width 
w(t) f3(D = 1) = 0.87 ± 0.04. The scaling behaviour of H (t) for larger systems with 
L = 262144 is plotted in fig. 3, while in fig. 4 the width w 2 (t) is shown for the same 
runs. The data were averaged over 40 independent runs. (This took 21.5 h CPU on 
one Cray Y-MP processor. The other data presented in this paper required about 
2000 h CPU on a IBM RS/6000 model 320 workstation.) The effective exponents 
(3{t) = log[w(2t)/w(t)]/ log2 (and analogously for H(t)) are shown in fig. 5a. A 
possibility to take into account corrections to scaling is to plot w 2 (2t) — w 2 (t) versus 
t. The corresponding effective exponents are shown in fig. 5b. The exponents 
determined by the width remain unchanged. Neither is there a clear improvement for 
the scaling of the height, although the plateau of f3(t) seems to be larger but also the 
fluctuations are bigger in fig. 5b. Note that the uncertainities in the last two points 
are much larger, so we do not expect that (3(t) increases to larger values. Thus, we 
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conclude (3(D = 1) = 0.88 ± 0.02. We have checked that this value is consistent with 
the case g = 2. The height-height correlation function C 2 (r, t) = ([h i+r (t) — hi(t)] 2 } 
is expected to scale in the same way as the width in eq.(3). Fig. 6 shows a scaling 
plot, where C(r, t) is divided by t 13 and r by t x l z . The best data collapse is achieved 
for ( ~ 1.25 and j3 ~ 0.88, in complete agreement with the above results. 

For large p = 0.95 and 0.97 the velocity is about v = 0.78 and 0.88, respectively 
(g = 2). From the interface width w(t) we find f3(p 3> p c ) = 0.25 ± 0.01, which is 
in agreement with the expected EW [14] regime far away from the depinning transi- 
tion. We were not able to estimate the exponent 9. Due to the large roughness the 
fluctuations of v become too large to fix the exponent 9 with a reasonable accuracy. 

b) 2+1 dimensions 

As in 1+1 dimensions we determine ( from the simulation of interfaces with 
different sizes until they get pinned. Fig. 7 shows the effective exponents C(-^) f° r 
(7 = 4, p = p c ~ 0.6416 and g = 6, p = p c ~ 0.74448 from which we estimate 
( = 0.75 ± 0.02. The scaling behaviour of the interface width for a system of size 
L 2 = 1024 2 is shown in fig.8 from which follows (3(D = 2) = 0.475 + 0.015. (The data 
were averaged over 15 independent runs.) The scaling of H(t) has to be considered 
with caution as will become clear from fig. 9 where the effective exponents f3(i) are 
plotted. For t > t x = AL Z , H(t) goes to a constant for p < p c or grows linearly 
for p > p c , while the scaling H(t) ~ t@ is expected for t a <^ t <^ t x where t a is the 
time where [3(t) achieves its asymptotic value. The plateau of [3(t) determined by 
H(t) for p > p c in fig. 9 (triangles) for 32 < t < 128 may be attributed to a very 
slow crossover because t a is of the same order as t x . Also for p = p c , (3(t) (circles) 
never achieves its asymptotic value. Thus, also the value for (3{t) determined from 
the scaling of the interface height H(t) in 1+1 dimensions has to be considered with 
caution. For g = 6 we find the same value f3 = 0.475 ± 0.015. (For g = 1 the velocity 
is nonzero even for p < q which would correspond to " negative driving forces" . This 
strange situation is perhaps due to the fact that the growth rule (5) does not allow 
hi to decrease. Thus we choose g such that p c > 1/2.) The best scaling plots of 
C(r,t) are achieved for the parameters ( = 0.74 and (3 = 0.47. ..0.48 consistent with 
the values given above. 
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In fig. 10 the velocity is plotted versus p — p c and with v ~ (p — p c ) e we get 
9 = 0.65 ± 0.05. Although the fluctuations of v are smaller in D = 2 than in D = 1 
the uncertainties remain rather big. For g = 4 and very small v, seems to be about 
0.68 but it then crosses over to 6 = 0.64. The exponent is very sensitive to the 
value of p c for small v. For larger v, on the other hand, it is not clear whether the 
scaling law v ~ (p — p c ) e still holds. It is known from the charge-density waves that 
the critical region is quite narrow [2]. So the error bar on 6 is larger than those on 
the roughness exponents ( and (3. 

IV. Discussion 

From the two roughness exponents ( and (3 we can calculate the dynamical expo- 
nent z = C//3: z(D = 1) = 1.42 ±0.03 and z(D = 2) = 1.58 ±0.04, i.e. the dynamics 
at the depinning transition is superdiffusive. The correlation length exponent v can 
be found from the general scaling relation v = 6/(z — () [15]. In 2+1 dimensions 
we get v{D = 2) = 0.8 ± 0.05. Simple arguments suggest [15] that there is a Harris 
criterion [21] for a sharp 2nd order depinning transition: \jv < (D + C)/2. In our 
case this relation is indeed fulfilled as an inequality. 

The simulations of the continuum equation (1) by Tang [18] yield a dynamical 
roughness exponent f3(D = 1) = 0.81 ± 0.02 which was determined from the scaling 
of the height H(t). This value is somewhat smaller than our result (3 = 0.88 ± 0.02. 
However, we have seen that it is rather difficult to settle the fluctuations of the 
width are large and the height shows slow crossover phenomena, especially for smaller 
systems. Another possibility is that due to the special discretization the dynamics 
of the lattice model is not described by the continuum equation (1), e.g. according 
to the growth rule eq.(5) there are only two velocities, zero and one, whereas in 
the continuum model the velocity depends continuously on the force. On the other 
hand the static roughness exponent £ ~ 1.25 agrees with the simulation [18] of 
eq.(l) but not with the Imry-Ma result [8] (£ = (4 - D)/Z for all D < 4) nor with 
the conjecture of Parisi [17] that £ = (4 - D)/2, (3 = (4 - L>)/4, 0=1. Our 
numerical values for the critical exponents support the perturbative renormalization 
group expansion in L> = 4 — e, e C 1 [15], because the extrapolation to e > 1 gives 
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a better estimate for e = 2 than for e = 3. It is interesting to note that the scaling 
relations z = (/ f3 = 2 — (e — £)/3 an d = 1 — |(e — C)/(2 — C) which were derived 
for small e [15] are perfectly fulfilled by our numerical results, although we have no 
arguments why they should be exact. By inserting our £ in the scaling relations for 
D = 1 we obtain 9(D = 1) ~ 0.22 and u(D = 1) ~ 1.33. 

Our roughness exponents are much larger than those of the EW model and the 
KPZ equation [22], i.e. quenched noise roughens an interface more than thermal 
noise. This difference corresponds to the observed crossover from the behaviour at 
the threshold {(3{D = 1) ~ 0.88) to that of large velocities (/3(D = 1) ~ 0.25) where 
the noise is short-range correlated in time. Next, we discuss the relation to the growth 
mechanism "pinning by directed percolation" [20], where the dominating noise is also 
quenched. There, if the height difference between two neighbouring columns exceeds 
a certain small value, the lower column grows, independent of random pinning forces. 
Thus, the difference between neighbouring columns is kept small and the growth is 
side-ways. It can be shown [18] that the model is in the same universality class as 
eq.(l) if a KPZ-term F(Vz) 2 /2 is added, which favours lateral growth. In the model 
presented here, however, only the discretized Laplacian, i.e. the sum of the height 
difference to all nearest neighbours enters the growth rule (5). Therefore, the motion 
is not side-ways, which corresponds to the absence of the gradient-squared term 
F(Vz) 2 /2 in eq.(l), and arbitrary large slopes are possible which can occur when 
the height difference to two neighbours has opposite sign. Indeed, our roughness 
exponents are larger than those of an interface pinned by a directed percolation cluster 
(CCD = 1) = (3{D = 1) ~ 0.63 at the threshold [20]). In 1+1 dimensions the exponent 
( ~ 1.25 > 1 is even unphysical. If one derives eq.(l) from a Hamiltonian [1,4] the 
gradient terms have to be assumed to be small. If however ( > I the neglection of 
overhangs and higher order gradients is no longer justified. It is therefore questionable 
to apply this simple model to real 1+1 dimensional systems. For the equilibrium case 
of the random field Ising model it is known [4] that d = D + 1 = 2 is the lower critical 
dimension di. At and below di no long range order can be established and the domain 
walls are convoluted in contradiction with the assumption of no overhangs. 

In 2+1 dimensions a recent simulation [5] of the motion of a domain wall in the 
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random field Ising model yields a roughness exponent ( = 0.67 ± 0.03 which was 
obtained by a plot of the height-height correlation function C(r) for systems up to 
linear size 300. In our measurements C(r) never reaches its asymtotic scaling with 
C ~ 0.75 due to finite time and finite size effects but gives effective exponents about 
2/3. It would be therefore very interesting to measure the width of a domain wall 
for different system sizes to check whether the domain walls in the random field Ising 
model of ref.[5] are correctly described by the model presented in this paper. 

V. Summary 

We have proposed a lattice model for the driven motion of an interface in a 
random medium. The critical exponents characterizing the depinning transition are 
determined in one and two interface dimensions. We believe that the model is in the 
same universality class as the continuum equation (1). Our results (e.g. ((D = 1) ~ 
1.25, C(D = 2) ~ 0.75) are inconsistent with the validity of an Imry-Ma argument 
[8] (C = (4 — -D)/3, D < 4) and with the approximation in eq.(l) to neglect the 
^-dependence of the random force r\ [17] (£ = (4 — D)/2), but they support the 
analytical results of the renormalization group for D = 4 — e, eC 1 [15] (( = e/3), 
because the deviations are smaller for e = 2 than for e = 3. It is of course desirable 
to investigate whether this difference decreases further when e = 1. 
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We determined the velocity exponent 9 in 1+1 dimemnsions by a simululation of 
large systems (up to L = 262144) to 6 = 0.25 ± 0.03. (This took 30 h CPU on one 
Cray Y-MP processor.) 
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Figure captions 

Fig.l The width w 2 (L) for g = 1, p = p c ~ 0.8004 (squares) and for g = 2, p = 
p c ~ 0.8748 (circles). The statistical uncertainities are smaller than the size of 
the symbols. 

Fig. 2 The effective roughness exponent ((L) for g = 1, p = p c (squares) and g = 2, p = 
0.8744 < p c (circles). For p < p c an interface is pinned at earlier times and one is 
able to simulate larger systems than for p = p c . 

Fig.3 Scaling of H(t) for L = 262144 and p = p c = 0.8004 (g = 1). 

Fig. 4 The width w 2 (t) for the same runs as in fig.3. 

Fig. 5 The effective exponents f3(t) determined from the height H(t) (circles) and from 
the width w 2 (t) (squares). In fig. 5b f3(t) is determined from a plot w 2 (2t) —w 2 (t) 
(and analogous H{2i) — H{t)), whereas in fig 5a the exponents are calculated from 
a simple plot w 2 (t) and H(t). 

Fig. 6 Best scaling plot of the height- height correlation function C(r, t), which is scaled 
by t/ 3 , (3 = 0.88 versus r ( = 1.25. The same plotting symbol is used for 

data at a given time t. In fig. 6a all data are plotted and in fig. 6b the crossover 
region is shown only. 

Fig. 7 The effective exponents £(L) for g = 4, p = p c ~ 0.6416 (squares) and g = 6, p = 
p c ~ 0.74448 (circles). 

Fig.8 The interface width w 2 (t) for L 2 = 1024 2 , g = 4, p = p c = 0.6416. 

Fig. 9 The effective exponents f3(t) determined from the plot of w 2 (t) (squares) and 
from H(t) (circles) for g = 4, p = p c = 0.6416. In addition, (3(t) from the height 
for the case p = 0.643 > p c (triangles) is plotted for comparison (see text). 

Fig. 10 Double logarithmic plot of the velocity versus p—p c for g = 4 (squares) and g = 6 
(circles). 
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